The Journal of Chemical Physics

The unfolding or melting temperature (TM) is a central quantity to characterize the stability of proteins and other biopolymers. The accurate prediction of protein melting temperatures by molecular mechanics force field simulations is highly desirable for many biophysical and biotechnological applications. Since the time scales for protein (un-)folding are hardly accessible in conventional MD (cMD) simulations, enhanced sampling techniques such as Temperature Replica Exchange Molecular Dynamics (TREMD) are typically employed. However, TREMD simulations are computationally very demanding especially if large temperature ranges need to be covered. Additionally, if the TM is initially unknown, setting up TREMD simulations is often challenging. To find the optimal initial conditions for such simulations, we describe their performance based on a theoretical model, which we validate on a minimalistic Markov Chain Monte Carlo (MCMC) simulation setup. In an effort to reduce the computational demand, we have investigated the possibility to use small sets of TREMD temperature ladders placed iteratively in the vicinity of a TM estimate. Different TREMD setups were extensively tested on the fast-folding protein Chignolin. We found that appropriate starting conformations lead to significantly faster convergence. Furthermore, we found that, in practice, combining multiple small temperature ladders can be advantageous in comparison to one single temperature ladder. Based on our findings, we formulate practical recommendations on how to set up TREMD for protein melting with optimal efficiency.

Molecular dynamics (MD) simulations are a powerful tool for investigating biomolecular dynamics underlying biological functions. However, the accessible spatiotemporal scales of conventional all-atom simulations remain limited by high computational costs. Coarse-graining reduces these costs by decreasing the number of interaction sites and enabling longer timesteps. In extreme cases, proteins are represented as single spherical particles; while such approximations facilitate cellular-scale simulations, they often sacrifice essential structural information, such as molecular shape and interaction anisotropy. Here, we present CGRig, a rigid-body protein model with residue-level interaction sites designed for long-time, large-scale simulations. In CGRig, each protein is treated as a single rigid-body embedding residue-level interaction sites. Its translational and rotational motions are described by the overdamped Langevin equation incorporating a shape-dependent friction matrix. Intermolecular interactions are calculated using G[o]-like native contact potentials, Debye-Huckel electrostatics, and volume exclusion. We validated that CGRig accurately reproduces the translational and rotational diffusion coefficients expected from the friction matrix for an isolated protein. For dimeric systems, the model successfully maintained native complex structures. Furthermore, two initially separated proteins converged into the correct complex with an association rate consistent with all-atom simulations. Notably, CGRig achieved a simulation performance exceeding 17 s/day for a 1,024-molecule system. These results demonstrate that CGRig provides an efficient framework for simulating protein assembly while retaining residue-level interaction specificity, making it a valuable tool for investigating large-scale biomolecular self-assembly.

Nuclear Pore Complexes (NPCs) are large protein complexes in eukaryotic cells that span the double-membrane of the nucleus and regulate bi-directional transport between nucleus and cytoplasm. T h e NPC core is lined by intrinsically disordered protein chains called nucleoporins (Nups) which form a selective barrier where large macromolecules (cargoes) need to bind to nuclear transport receptors (NTRs) such as Karyopherins (Kaps) to cross. Previous experimental results have suggested that not only Nups but Kaps, too, are important in the transport process of other NTRs/NTR-cargo complexes. In this work, we assess the role of Kaps in the transport of other NTRs (specifically, NTF2s) through the NPC, a process referred to as the "Kap-centric transport model". Here, using coarse-grained MD simulation we show that Kaps are able to direct NTF2s into the Nup meshwork, which leads to their increased flow. Our results also suggest that NTRs follow specific lanes inside the pore to maximize efficient transport.

There is a growing demand for molecular dynamics simulations to explore longer timescale behavior of giant protein-DNA complexes such as chromatin. To address this need, we extended OpenCafeMol, a GPU-accelerated residue-level coarse-grained molecular dynamics simulator originally developed for proteins and lipids, to support 3SPN.2 and 3SPN.2C DNA models. We also implemented a hydrogen-bond-type many-body potential to model DNA-protein interactions more accurately. To further improve computational efficiency, we introduced a localized scheme for calculating base-pairing and cross-stacking interactions. Benchmark tests show that OpenCafeMol on a single GPU achieves up to 200-fold speed-up for DNA-only systems and up to 100-fold speed-up for DNA-protein complexes compared to CPU-based simulations. To demonstrate the capability of our implementation for long-timescale biological processes, we simulated an archaeal SMC-ScpA complex undergoing DNA translocation via segment capture (a proposed mechanism for DNA loop extrusion) in the presence of a DNA-bound obstacle. We observed continuous captured-loop growth accompanied by obstacle bypass within the segment capture framework.

Accurate ion force fields are essential for molecular dynamics simulations of biomolecular systems, particularly in combination with modern water models such as OPC. While OPC water improves the description of bulk water and biomolecules, the transferability of existing ion force fields to this model remains an open question. Here, we systematically assess the transferability of monovalent and divalent ion force field parameters (Li+, Na+, K+, Cs+, Mg2+,Ca2+, Sr2+, Ba2+, Cl- and Br-) to OPC water by comparing single-ion and ion-pairing properties with experimental data. Our analysis reveals that no single literature parameter set provides accurate results for all ions when directly transferred to OPC water. We hence introduce the MS/G-LB(OPC) force field, which combines Mamatkulov-Schwierz-Grotz cation parameters with Loche-Bonthuis anion parameters. MS/G-LB(OPC) reproduces hydration free energies, first-shell structural properties and activity derivatives at low salt concentrations. Our results demonstrate that transferring ion parameters to OPC can lead to significant and ion-specific deviations from experimental data, making careful validation essential. At the same time, the systematic transfer and combination of ion parameters from existing force fields can provide a practical and computationally efficient alternative to full reparameterization. MS/G-LB(OPC) is available at https://git.rz.uni-augsburg.de/cbio-gitpub/opc-ion-force-fields.

The synaptic cleft between neighboring neurons is the site of neurotransmitter-mediated communication that underlies normal brain function, including learning and memory. When an action potential reaches the presynaptic terminal, released neurotransmitters cross the cleft under the combined influence of diffusion and electrical forces to activate postsynaptic receptors. Despite this, synaptic-cleft transport is commonly modeled using a pure diffusion model, neglecting electrical drift. Here, we quantify the relative contributions of diffusion and electrical terms in the Poisson-Nernst-Planck (PNP) framework and assess whether the pure diffusion approximation is adequate. We solve the full PNP system in a three-dimensional computational model of the synaptic cleft at nanometer-scale resolution, tracking five ionic species (Na+, K+, Ca2+, Cl-, Glu-) with full spatial and temporal detail. Solutions are compared directly with those of the pure diffusion (D) model. The D and PNP models produce markedly different ionic concentration fields. Analysis of ionic fluxes confirms that diffusive and electrical contributions are of comparable magnitude across all species. These discrepancies are robust across parameter variations, including the number of AMPA receptors, the amount of released glutamate, the cleft height, and the cleft diffusion coefficient, and are amplified as the number of AMPA receptors increases, the cleft becomes narrower or diffusion more restricted. The quantitative and qualitative differences between the pure D model and the full PNP model demonstrate that neglecting electrical forces in the synaptic cleft has consequences. These discrepancies are large enough to alter the predicted dynamics and biological interpretation of synaptic transmission, establishing that accurate computation of ionic concentrations in the synaptic cleft requires the full PNP equations.

Targeted protein degradation is an emerging approach that utilizes cellular degradation pathways to inhibit a target protein. Small molecules such as molecular glues or PROTACs can be used to mediate the formation of a ternary complex with an E3 ligase and the target protein, which can dramatically enhance the degradation process. This approach is promising for cancer therapy, where degradation of oncogenic proteins can lead to cancer cell toxicity. To design new molecular glues, it is important to develop methods that predict how well a given molecule stabilizes a protein-protein interaction. However, conventional molecular dynamics simulations face challenges in capturing the long-timescale binding and unbinding events that would be used to evaluate this stabilization. In this study, we developed a strategy that allows us to evaluate the stability of protein-protein interactions in the presence of a glue molecule using weighted ensemble simulations in combination with weakened protein-protein interactions. Using this strategy, we generated unbinding trajectories of the DCAF15-RBM39 system with small molecules E7820, Indisulam, and several other Indisulam analogs. We were able to observe distinctly different behaviors between systems with different glues, which was in agreement with their reported EC50 values. We believe this approach could aid drug discovery efforts by expanding the set of druggable targets and improving the success rate of molecular glue development.

The flexibility of DNA and RNA is known to play a central role in numerous biological processes, including chromatin organization and gene regulation. While a wide range of computational approaches have been developed to investigate the conformational dynamics and flexibility of proteins, analogous methods for nucleic acids remain comparatively underexplored. Elastic Network Models (ENMs) - coarse-grained mechanical representations in which macromolecules are modeled as networks of nodes connected by elastic springs - have been successfully applied to proteins, often allowing to capture experimentally observed conformational changes through a small number of harmonic normal modes. Building on a previously validated three-bead ENM for RNA, here we introduce edENM, an essential dynamics-refined ENM for DNA, RNA, and protein-nucleic acid complexes, parametrized using a diverse set of Molecular Dynamics simulations. The vibrational modes of the new edENM show good agreement with NMR data and experimental ensembles, while avoiding the unrealistic and localized deformability of previous ENM parametrizations. Additionally, we integrated this new edENM into eBDIMS, a Brownian Dynamics-based framework that enables the simulation of large-scale and anharmonic conformational transitions in protein assemblies. In this way, we are now able to explore functional motions in large protein-nucleic acid complexes such as chromatin subunits and ribosomes.

Reconstructing all-atom (AA) structures from highly coarse-grained (HCG) models remains a significant challenge in multiscale molecular dynamics (MD) simulations, particularly for mesoscale biomolecular assemblies that are beyond the reach of conventional MD methods. Building upon ProNet Backmapping, a neural-network-based thermodynamically consistent approach, we introduce a progressive backmapping framework that reconstructs AA models in a stepwise manner across neighboring resolutions, for example, from a 3-residue-per-site HCG model to a 1-residue-per-site model, then to an AA model. This progressive backmapping method achieves high accuracy across a wide range of proteins and effectively reconstructs flexible linkers in multidomain architectures. Moreover, it supports hierarchical reconstruction of complex protein assemblies, including multiple virus-like particles spanning tens of nanometers and containing hundreds of subunits. Using this framework, we demonstrate--for the first time--the ability to hierarchically backmap entire viral assemblies from HCG to full AA resolution, covering at least three different resolutions. Overall, our method provides a scalable framework for incorporating atomistic detail into mesoscale simulations of complex systems across many applications in chemistry and biology. Table of contents figure O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=70 SRC="FIGDIR/small/709104v1_ufig1.gif" ALT="Figure 1"> View larger version (45K): org.highwire.dtl.DTLVardef@4af423org.highwire.dtl.DTLVardef@e2669borg.highwire.dtl.DTLVardef@1be80eforg.highwire.dtl.DTLVardef@2e679_HPS_FORMAT_FIGEXP M_FIG C_FIG

Accessing deeply buried binding sites remains a major challenge in structure-based drug discovery, where accurate description of both protein dynamics and ligand binding pathways is required. Funnel metadynamics enables simulation of complete binding processes but is computationally demanding at the all-atom resolution. By adopting the Martini 3 force field, coarse-grained funnel metadynamics (CG-FMD) substantially reduces computational requirements while retaining enhanced sampling capabilities. In this work, we assess the capability of CG-FMD to model ligand recognition at the deeply buried colchicinoids site of the tubulin {beta}-heterodimer, a multisite protein of strategic importance. We investigated the binding of colchicine, podophyllotoxin and combretastatin-A4, recovering free energy profiles with improved statistical convergence compared to AA-FMD and comparable to experimental references. In particular CG-FMD binding free energies present mean absolute errors between 3 and 10 kJ mol-1. These results propose CG-FMD as an efficient, physics-based framework for probing ligand binding to challenging sites.

Cellular adaptation is inherently nonstationary processes with complex stochastic dynamics1-5. Despite remarkable progress in quantitative biology6-11, a quantitative understanding of the cell adaptation dynamics in terms of the underlying cellular network remains elusive. Here, we present the next-generation chemical dynamics model and theory for cellular networks, providing an effective, quantitative description of the adaptive gene expression dynamics in living cells responding to external stimuli. Unlike conventional kinetics, chemical dynamics of cellular network modules are characterized by their reaction-time distributions, rather than by rate coefficients12. For a general model of cell signal transduction and adaptive gene expression, we derive exact analytical expressions for the time-dependent mean and variance of protein numbers produced in response to external stimuli, validated by accurate stochastic simulations. These results provide a unified, quantitative explanation of the stochastic responses of diverse E. coli genes to antibiotic stress and transcriptional induction. Our analysis reveals existence of a general quadratic relationship between the mean and variance of activation times across diverse genes. The gene activation process influences transient dynamics of downstream protein levels, but not their steady-state levels. In contrast, post-translational maturation process affects both transient dynamics and steady-state variability of mature protein levels. This finding indicates that the gene expression variability measured by fluorescent reporter proteins depends on the maturation time of the reporters. This work suggests a new direction for the development of digital twins of living cells.

Biomolecular condensates formed by complex coacervation of highly charged proteins provide a powerful framework to understand how microscopic interactions give rise to macroscopic material properties. Atomistic molecular dynamics simulations provide detailed insights but remain limited in accesing the spatio-temporal scales relevant for condensate behavior. Here, we use the residue-level coarse-grained Mpipi-Recharged model to investigate condensates formed by ProT and positively charged partners, including histone H1, protamine, poly-lysine, and poly-arginine. Material properties, in this context, provide a stringent experimental benchamark for coarse-grained models. Our model reproduces salt-dependent phase behavior, protein binding affinities, and sequence-specific stability trends in agreement with in vitro experiments, despite the fact that material properties were not included in the model parametrization. We then establish a direct link between protein dynamics and macroscopic material properties by quantifying monomeric diffusion, conformational reconfiguration, and translational mobility within the dense phase, and relating these to condensate viscosity. By comparing dynamics across dense and dilute phases, we uncover a pronounced length scale-dependent behavior. While residue-level binding and unbinding events remain equally fast in both phases, protein reconfiguration time and self-diffusion are significantly slowed down within the condensates. This decoupling reveals how fast intermolecular interactions coexist with slow mesoscale condensate dynamics depending on the molecular length scale. Together, our results establish a predictive framework that links encoded sequence intermolecular forces and multiscale dynamics to the emergent material properties of complex biomolecular condensates.

The interactions of droplets and filaments can lead to mutual deformations and complex combined behavior. Such interactions also occur within the cell, where biomolecular condensates, distinct liquid phases often composed of proteins, have been observed to structure and affect the organization of the cytoskeleton. In particular, biomolecular condensates have been shown to undergo characteristic deformations when cytoskeletal filaments are fully embedded within them. However, a full understanding of the underlying physical mechanisms is still missing. Here, we combine experiments with coarse-grained molecular dynamics simulations and analytical models to uncover the physical mechanisms that define emerging shapes of droplets containing filaments. We find that the surface tension of the liquid phase and the bending energy of the filament(s) suffice to accurately capture emerging shapes if the length of the filament is small compared to the liquid volume. As the volume fraction of filament(s) increases, wetting effects become increasingly important, setting physical constraints within which surface and bending energies compete to define the droplet shapes. We find that mutual deformations of condensate and filament extend accessible shapes beyond classical stability considerations, leading to structuring and entrapment of contained filaments. Shape deformations may further affect ripening dynamics that favor certain geometries. Our findings provide a physical framework for a better understanding of the possible roles of biomolecular condensates in cytoskeletal organization.

We present a fast, unconditionally energy-stable numerical scheme for simulating vesicle deformation under osmotic pressure using a phase-field approach. The model couples an Allen-Cahn equation for the biomembrane interface with a variable-mobility Cahn-Hilliard equation governing mass exchange across the membrane. Classical approaches, including nonlinear multigrid and Multiple Scalar Auxiliary Variable (MSAV) methods, require iterative solution of variable-coefficient systems at each time step, resulting in substantial computational cost. We introduce a constant-coefficient MSAV (CC-MSAV) scheme that incorporates stabilization directly into the Cahn-Hilliard evolution equation rather than the chemical potential. This reformulation yields fully decoupled constant-coefficient elliptic problems solvable via fast discrete cosine transform (DCT), eliminating iterative solvers entirely. The method achieves O(N2 log N) complexity per time step while preserving unconditional energy stability and discrete mass conservation. Numerical experiments verify second-order temporal and spatial accuracy, mass conservation to relative errors below 5 x 10-11, and close agreement with nonlinear multigrid benchmarks. On grids with N [≥] 2048, CC-MSAV achieves 6-15x overall speedup compared to classical MSAV with optimized preconditioning, while the dominant Cahn-Hilliard subsystem is accelerated by up to two orders of magnitude. These efficiency gains, achieved without sacrificing accuracy, make CC-MSAV particularly well-suited for large-scale simulations of vesicle dynamics.

We develop a self-consistent free-energy framework in which membrane shape and osmotic pressure are determined simultaneously in a finite reservoir by minimizing bending elasticity and solute entropy. Solute conservation makes osmotic pressure a thermodynamic variable rather than an externally prescribed parameter, producing a nonlinear coupling between membrane mechanics and solvent entropy. This coupling modifies the classical stability condition for spherical vesicles: instability emerges from global free-energy competition rather than the linear Helfrich stability criterion. The resulting critical pressures differ by orders of magnitude from Helfrich predictions and agree with simulations for small and large unilamellar vesicles. The framework is relevant to cellular environments involving biomolecular condensate confinement as well as synthetic vesicles and the development of osmotic-pressure-driven encapsulation platforms.

Disordered proteins can form biomolecular condensates by demixing from their environment, enabling reversible compartmentalisation of cellular components in the form of membraneless organelles. Multivalent interactions are essential for this type of phase separation behaviour, and for disordered proteins, the potential for multivalent interactions is encoded in the sequence composition and patterning. Mutational studies have been instrumental in helping elucidate this sequence grammar by perturbing the amino acid sequence and quantifying the resulting changes in the driving force for phase separation. While such studies have provided a detailed and predictive understanding of the driving forces for phase separation, they strictly do not inform on the nature of the interactions that drive phase separation. Here, we propose using double mutant cycles to explore molecular interactions and their contributions to condensate properties more directly. We explore the applicability of double mutant cycles for different types of interactions in condensates formed by the low-complexity domain of hnRNPA1 using coarse-grained molecular dynamics simulations. We find that the interactions between arginine and tyrosine residues, as well as between aromatic residues, contribute mostly additively to the propensity for phase separation. However, for the interactions between charged residues, we find that--in an interplay with the net charge of the protein--there is a measurable non-additive contribution to the phase separation propensity. Based on our results, we envisage that double mutant cycles could provide additional insights into protein phase separation, thus expanding the understanding of the sequence grammar and the underlying molecular interactions.

The spatial architecture and mechanical rigidity of polysomes are crucial determinants of translational efficiency and mRNA stability. In this study, we investigate the conformational statistics of an mRNA backbone decorated with high-density ribosomes at varying densities using large-scale, extensive molecular dynamics simulations based on the Kremer-Grest bead-spring model. To address the extreme spatial asymmetry between mRNA monomers and ribosomes, we used an efficient tree-based neighbour list algorithm, enabling the analysis of mRNA chains up to N = 4, 969. Our results demonstrate that the excluded volume of massive ribosomes induces a significant and robust expansion of the scaling exponent v from 0.59 to approximately 0.7. In the conformation of mRNA, this shift translates to a self-induced dimensional reduction from a three-dimensional random coil toward a stretched, a quasi-two-dimensional architecture at biologically relevant scales. Such a transition is further evidenced by a periodic "regain" of the bond-bond correlation function C(n) at ribosome attachment sites, indicating a geometric alignment absent in standard homopolymers. These findings reveal that the geometric crowding of ribosomes itself provides a robust physical prerequisite for the formation of higher-order polysome architectures, bridging the gap between polymer physics and structural properties of mRNA during translation.

Coarse-grained (CG) modeling of DNA has been gathering steam in recent years due to the limitations imposed by current computer hardware, and deficiencies of atomistic force-fields. Specifically, CG simulations have emerged as a potent tool for exploring complex energy landscapes underlying biochemical processes, such as DNA transcription, as well as material design based on programmable self-assembly. In this chapter, we illustrate how the Three Interaction Site (TIS) model for DNA, a robust coarse-graining framework, can be used to study the folding landscape of DNA hairpins. We show that despite its simplicity, the TIS-DNA model quantitatively describes the hairpin folding thermodynamics and recapitulates many features of the kinetics, including the multiplicity of pathways. The free energy landscape exhibits single-funnel character with a distinct bias towards the folded state. It is likely that folding initiates through non-specific collapse of the DNA chain, involving multiple excursions on the energy landscape, until the opposing strands are approximately aligned. Subsequently, the loop region becomes more ordered, and after the first native-contact nucleates, the rest of the process becomes essentially downhill.

Metabolites are ubiquitous in all living cells and are essential mediators of biochemical processes, serving either as substrates or as cofactors to enable the reactions. Capturing this diversity in computational workflows is important for allowing realistic simulations of the cytoplasm. Coarse-grained molecular dynamics enables the simulation of large scale systems up to the level of whole-cells, but is limited by the availability of refined parameters for all possible components in the system. In this work, we describe the parameterization of 186 common metabolites found in bacteria and eukaryotes within the framework of the Martini 3 force field. To showcase the behavior of Martini metabolites in a biological setting, we report simulations of protein-ligand binding and membrane permeation. The establishment of a Martini metabolome enables high-throughput simulations of metabolites interacting with other biomolecules, and opens the way for simulations of realistic cellular environments.

Molecular dynamics (MD) simulations have yielded important insights into ion conduction in potassium channels, but quantitative comparison with electrophysiological experiments remains challenging. Due to their high computational cost, MD simulations are typically performed at membrane potentials well above physiological values, and at only a limited number of voltages. Since current-voltage relationships are not necessarily linear, this limits direct comparison between simulations and experiments. Here, we introduce a method to estimate the current-voltage characteristics of ion channels from Markov state models (MSMs) constructed from MD simulations performed at only a few membrane potentials. Time-discrete MSMs of ion conduction are converted into continuous-time rate matrices, whose voltage dependence is modelled using a rate theory formulation with free energy barriers depending on membrane potential. This approach enables the prediction of channel currents over a wide voltage range without additional simulations. We validated the method using MD simulations of the potassium channels KcsA and MthK. In both cases, the currents predicted at low membrane potentials are in good agreement with those obtained directly from MD simulations, demonstrating the robustness and efficiency of the approach.